# coding: utf-8
import io
import sys
import numpy as np
import math
from pathlib import Path
import heapq
import cv2
import PIL
import scipy.signal as signal
import matplotlib.pyplot as plt
import mmcv

__all__ = ['calculate_cobb', 'calculate_line_angle', 'calculate_kps_angles']


def calculate_line_angle(kp:np.ndarray):
    k = (kp[1][1] - kp[0][1]) / (kp[1][0] - kp[0][0])
    angle = np.arctan(k) / np.pi * 180
    return angle


def calculate_kps_angles(kps:np.ndarray):
    """
    计算终板与水平线的夹角
    :param kps:
    1,2
    3,4
    5,6
    :return:
    """
    # 每块骨骼左右的中心点
    angles = []
    # for j in range(int(len(kps) / 4)):
    #     middle_left = ((kps[j * 4 + 2] + kps[j * 4]) / 2).astype(int)
    #     middle_right = ((kps[j * 4 + 3] + kps[j * 4 + 1]) / 2).astype(int)
    #     angles.append(calculate_line_angle(np.vstack((middle_left, middle_right))))

    for j in range(int(len(kps) / 2)):
        middle_left = (kps[j * 2]).astype(int)
        middle_right = (kps[j * 2 + 1]).astype(int)
        angles.append(calculate_line_angle(np.vstack((middle_left, middle_right))))

    return angles


def plot_angles_bar(angles, peaks, valleys):
    num = len(angles)
    x = np.arange(num)
    vertebra_x = ['T1', 'T1','T2', 'T2','T3', 'T3','T4','T4', 'T5', 'T5','T6', 'T6',
                  'T7','T7', 'T8', 'T8','T9', 'T9','T10', 'T10','T11', 'T11','T12','T12',
                  'L1','L1','L2','L2', 'L3','L3', 'L4', 'L4','L5', 'L5']
    # plt.bar(x, height=angles, width=1, color="green")
    plt.figure(figsize=(10, 10))

    # plt.rcParams['savefig.dpi'] = 600  # 图片像素
    # plt.rcParams['figure.dpi'] = 300  # 分辨率

    plt.plot(x, angles, linestyle='dashed',linewidth=0.5,color='green',marker='.',label='angles line')
    plt.xlabel("index")
    plt.ylabel("angle")
    plt.title("angles bar")
    plt.xlim((0, num))
    plt.xticks(x, vertebra_x)

    plt.hlines(0, 0, 33, color="blue", linestyle='--')
    # plt.draw()
    # plt.pause(10)
    plt.scatter(np.array(peaks)[:, 0], np.array(peaks)[:, 1], color='blue')
    plt.scatter(np.array(valleys)[:, 0], np.array(valleys)[:, 1], color='red')

    buffer = io.BytesIO()
    plt.savefig(buffer, format='png')
    buffer.seek(0)
    image = PIL.Image.open(buffer)

    ar = np.asarray(image)
    # cv2.imshow("demo", ar)
    # cv2.waitKey(0)
    buffer.close()
    plt.close()
    return ar


def find_peak_velley(angles, delta=0.3, x = None):
    # angles = np.array(angles)
    # peaks, _ = signal.find_peaks(angles, distance=3)
    # peaks = peaks.astype(int)
    # plt.plot(angles)
    # plt.plot(peaks, angles[peaks], "x")
    # plt.show()
    """
    Converted from MATLAB script at http://billauer.co.il/peakdet.html

    Returns two arrays

    function [maxtab, mintab]=peakdet(v, delta, x)
    %PEAKDET Detect peaks in a vector
    %        [MAXTAB, MINTAB] = PEAKDET(V, DELTA) finds the local
    %        maxima and minima ("peaks") in the vector V.
    %        MAXTAB and MINTAB consists of two columns. Column 1
    %        contains indices in V, and column 2 the found values.
    %
    %        With [MAXTAB, MINTAB] = PEAKDET(V, DELTA, X) the indices
    %        in MAXTAB and MINTAB are replaced with the corresponding
    %        X-values.
    %
    %        A point is considered a maximum peak if it has the maximal
    %        value, and was preceded (to the left) by a value lower by
    %        DELTA.

    % Eli Billauer, 3.4.05 (Explicitly not copyrighted).
    % This function is released to the public domain; Any use is allowed.

    """
    maxtab = []
    mintab = []

    if x is None:
        x = np.arange(len(angles))

    v = np.asarray(angles)

    if len(v) != len(x):
        sys.exit('Input vectors v and x must have same length')

    if not np.isscalar(delta):
        sys.exit('Input argument delta must be a scalar')

    if delta <= 0:
        sys.exit('Input argument delta must be positive')

    mn, mx = np.Inf, -np.Inf
    mnpos, mxpos = np.NaN, np.NaN

    lookformax = True

    for i in np.arange(len(v)):
        this = v[i]
        if this > mx:
            mx = this
            mxpos = x[i]
        if this < mn:
            mn = this
            mnpos = x[i]

        if lookformax:
            if this < mx - delta:
                maxtab.append((mxpos, mx))
                mn = this
                mnpos = x[i]
                lookformax = False
        else:
            if this > mn + delta:
                mintab.append((mnpos, mn))
                mx = this
                mxpos = x[i]
                lookformax = True

    return np.array(maxtab), np.array(mintab)
    # plt.plot(angles)
    # plt.scatter(np.array(maxtab)[:, 0], np.array(maxtab)[:, 1], color='blue')
    # plt.scatter(np.array(mintab)[:, 0], np.array(mintab)[:, 1], color='red')
    # plt.show()


def draw_peak_width(angles):
    angles = np.array(angles)
    peaks, _ = signal.find_peaks(angles)
    results_half = signal.peak_widths(angles, peaks, rel_height=1)
    plt.plot(angles)
    plt.hlines(*results_half[1:], color="C2")

    # neg_angles = angles * -1
    # peaks, _ = signal.find_peaks(neg_angles)
    # results_half = signal.peak_widths(neg_angles, peaks, rel_height=1)
    # #plt.plot(neg_angles)
    # plt.hlines(*results_half[1:], color="C2")
    plt.show()


def less_than(a, b):
    if a < b:
        return True
    else:
        return False


def great_than(a, b):
    if a > b:
        return True
    else:
        return False


def filter_peak_velley(angles:np.ndarray, peaks:np.ndarray, valleys:np.ndarray):
    cache_peaks = np.copy(peaks)
    cache_valleys = np.copy(valleys)
    num = len(angles)
    argmax_peak = np.argmax(cache_peaks[:, 1])
    max_peak_index = cache_peaks[argmax_peak][0]
    argmin_valley = np.argmin(cache_valleys[:,1])
    min_valley_index = cache_valleys[argmin_valley][0]
    filter_peaks_indexs = []
    filter_valleys_indexs = []
    filter_peaks_indexs.append(int(max_peak_index))
    filter_valleys_indexs.append(int(min_valley_index))

    def search_peak_valley_func(search_list,
                                cache_list,
                                filter_indexs,
                                offset=0,
                                argc_func=np.argmin):
        cur_min_valley_index = 1
        if len(search_list) != 0:
            cur_argmin_valley = argc_func(search_list[:, 1])
            cur_min_valley_index = int(search_list[cur_argmin_valley][0])
            if len(cache_list) != 0:
                left_min_valley_index = argc_func(cache_list) + offset
                if argc_func == np.argmin:
                    compare_func = less_than
                else:
                    compare_func = great_than
                if compare_func(angles[left_min_valley_index], search_list[cur_argmin_valley][1]):
                    cur_min_valley_index = left_min_valley_index
            filter_indexs.append(cur_min_valley_index)
        else:
            if len(cache_list) != 0:
                left_min_valley_index = argc_func(cache_list) + offset
                cur_min_valley_index = left_min_valley_index
                filter_indexs.append(cur_min_valley_index)
        return cur_min_valley_index

    # 1.从最大波峰和最大波谷 往两侧寻找代偿弯
    # 2.从两边同时找1个波峰和1个波谷,放入波峰、波谷缓存
    # 3.对缓存,进行过滤，取差值最大的三个配对值.

    # 判断最大波峰、最大波谷顺序
    if max_peak_index < min_valley_index:
        # 波峰在波谷前面
        # 1.寻找波峰前面
        # 1.1寻找波峰前面第一个波谷
        search_valleys = cache_valleys[np.where(cache_valleys[::, 0] < max_peak_index)[0], :]
        left_cache = angles[0:int(max_peak_index)]
        cur_min_valley_index = search_peak_valley_func(search_valleys,
                                                       left_cache,
                                                       filter_valleys_indexs,
                                                       0,
                                                       np.argmin)

        # 1.2寻找紧挨着后面的第一个波峰
        search_peaks = cache_peaks[np.where(cache_peaks[::, 0] < cur_min_valley_index)[0], :]
        left_cache = angles[0: int(cur_min_valley_index)]
        cur_max_peak_index = search_peak_valley_func(search_peaks,
                                                     left_cache,
                                                     filter_peaks_indexs,
                                                     0,
                                                     np.argmax)

        # 2.寻找波谷后面
        # 2.1 寻找波谷后面第一个波峰
        search_peaks = cache_peaks[np.where(cache_peaks[::, 0] > min_valley_index)[0], :]
        right_cache = angles[int(min_valley_index+1)::]
        cur_max_peak_index = search_peak_valley_func(search_peaks,
                                                     right_cache,
                                                     filter_peaks_indexs,
                                                     int(min_valley_index + 1),
                                                     np.argmax)

        # 2.2 寻找紧挨着后面第一个波谷
        search_valleys = cache_valleys[np.where(cache_valleys[::, 0] > cur_max_peak_index)[0], :]
        right_cache = angles[(cur_max_peak_index+1)::]
        cur_min_valley_index = search_peak_valley_func(search_valleys,
                                                       right_cache,
                                                       filter_valleys_indexs,
                                                       int(cur_max_peak_index + 1),
                                                       np.argmin)
    else:
        # 波峰在波谷后面
        # 1.寻找波谷前面
        # 1.1 寻找波谷前面第一个波峰
        search_peaks = cache_peaks[np.where(cache_peaks[::, 0] < min_valley_index)[0], :]
        left_cache = angles[0: int(min_valley_index)]
        cur_max_peak_index = search_peak_valley_func(search_peaks,
                                                     left_cache,
                                                     filter_peaks_indexs,
                                                     0,
                                                     np.argmax)

        # 1.2 寻找紧挨着前面第一个波谷
        search_valleys = cache_valleys[np.where(cache_valleys[::, 0] < cur_max_peak_index)[0], :]
        left_cache = angles[0: cur_max_peak_index]
        cur_min_valley_index = search_peak_valley_func(search_valleys,
                                                       left_cache,
                                                       filter_valleys_indexs,
                                                       0,
                                                       np.argmin)

        # 2.1寻找波峰后面第一个波谷
        search_valleys = cache_valleys[np.where(cache_valleys[::, 0] > max_peak_index)[0], :]
        right_cache = angles[int(max_peak_index+1)::]
        cur_min_valley_index = search_peak_valley_func(search_valleys,
                                                       right_cache,
                                                       filter_valleys_indexs,
                                                       int(max_peak_index + 1),
                                                       np.argmin)
        # # 2.2寻找紧挨着后面的第一个波峰
        search_peaks = cache_peaks[np.where(cache_peaks[::, 0] > cur_min_valley_index)[0], :]
        right_cache = angles[int(cur_min_valley_index + 1)::]
        cur_max_peak_index = search_peak_valley_func(search_peaks,
                                                     right_cache,
                                                     filter_peaks_indexs,
                                                     int(cur_min_valley_index + 1),
                                                     np.argmax)

    filter_peaks = []
    filter_valleys = []
    for v in filter_peaks_indexs:
        filter_peaks.append([v, angles[v]])
    for v in filter_valleys_indexs:
        filter_valleys.append([v, angles[v]])
    filter_peaks = np.array(filter_peaks)
    filter_valleys = np.array(filter_valleys)

    # 过滤：
    # 1.**两侧波峰波谷间距最小的排除**
    # 2.两侧排除没有跨越0的波峰波谷对
    # 最后保证的顶点《=4
    sort_vertexes: np.ndarray = np.sort(np.hstack((filter_peaks_indexs, filter_valleys_indexs)))
    diff_vertex = []
    for i in range(len(sort_vertexes) - 1):
        diff_vertex.append(abs(angles[sort_vertexes[i]] - angles[sort_vertexes[i+1]]))

    top_index = list(map(diff_vertex.index, heapq.nlargest(3, diff_vertex)))
    top_index.sort()
    top_index.extend([top_index[-1]+1])
    top_vertexes = sort_vertexes[top_index]

    finaly_peaks = set(top_vertexes.tolist()) & set(np.array(filter_peaks)[:, 0].tolist())
    finaly_valleys = set(top_vertexes.tolist()) & set(np.array(filter_valleys)[:, 0].tolist())
    finaly_peaks = np.array(list(finaly_peaks)).astype(filter_peaks.dtype)
    finaly_valleys = np.array(list(finaly_valleys)).astype(filter_valleys.dtype)

    return filter_peaks[np.in1d(filter_peaks[:, 0], finaly_peaks)], \
           filter_valleys[np.in1d(filter_valleys[:, 0], finaly_valleys)]


def max_bend(sort_vertexes):
    extre_list = sort_vertexes[:, 1].tolist()
    diff_list = []
    for i in range(len(extre_list) - 1):
        diff_list.append(abs(extre_list[i + 1] - extre_list[i]))
    # 解决某些情况下最大波峰或者波谷不会生成主弯的情况
    top_index = list(map(diff_list.index, heapq.nlargest(1, diff_list)))[0]
    extre_peak_valley = sort_vertexes[top_index:(top_index + 2), :]
    return extre_peak_valley, top_index


def calculate_cobb_by_angle(pts):
    pts = np.asarray(pts, np.float32)  # 68 x 2
    num_pts = pts.shape[0]  # number of points, 68
    # vnum = num_pts // 4 - 1
    # mid_p_v = (pts[0::2, :] + pts[1::2, :]) / 2  # 34 x 2
    # mid_p = []
    # for i in range(0, num_pts, 4):
    #     pt1 = (pts[i, :] + pts[i + 2, :]) / 2
    #     pt2 = (pts[i + 1, :] + pts[i + 3, :]) / 2
    #     mid_p.append(pt1)
    #     mid_p.append(pt2)
    # mid_p = np.asarray(mid_p, np.float32)  # 34 x 2
    # vec_m = mid_p[1::2, :] - mid_p[0::2, :]  # 17 x 2
    vec_m = pts[1::2, :] - pts[0::2, :]     # 左右点的连线
    dot_v = np.matmul(vec_m, np.transpose(vec_m))  # 17 x 17
    mod_v = np.sqrt(np.sum(vec_m ** 2, axis=1))[:, np.newaxis]  # 17 x 1
    mod_v = np.matmul(mod_v, np.transpose(mod_v))  # 17 x 17
    cosine_angles = np.clip(dot_v / mod_v, a_min=0., a_max=1.)
    angles = np.arccos(cosine_angles)  # 17 x 17
    pos1 = np.argmax(angles, axis=1)
    maxt = np.amax(angles, axis=1)
    pos2 = np.argmax(maxt)
    cobb_angle1 = np.amax(maxt)
    cobb_angle1 = cobb_angle1 / np.pi * 180

    angle2 = angles[pos2, :(pos2 + 1)]
    cobb_angle2 = np.max(angle2)
    pos1_1 = np.argmax(angle2)
    cobb_angle2 = cobb_angle2 / np.pi * 180

    return [pos2, pos1[pos2], pos1_1]


def calculate_cobb(kps:np.ndarray):
    """
    根据脊椎的68个关键点计算PT MT TL
    :param kps:keypoints shape [68, 2]
    :return:

    bend_angles: {'pt': {'angle':30, 'line_index': [5, 10]}, ....}

    """
    thoracic_index = {'PT':(0, 5), 'MT':(6, 12), 'TL/L':(13, 17)}
    thoracic_keys = ['PT', 'MT', 'TL/L']
    # 1.计算夹角
    angles = calculate_kps_angles(kps)
    # 2.找出所有的波峰波谷
    peaks, valleys = find_peak_velley(angles, 0.3)
    # 3.过滤波峰波谷
    peaks, valleys = filter_peak_velley(angles, peaks, valleys)
    sort_vertexes: np.ndarray = np.vstack((peaks, valleys))
    sort_vertexes = sort_vertexes[sort_vertexes[:, 0].argsort()]

    # 二次校验
    assert sort_vertexes.shape[0] == 4
    if int(sort_vertexes[1, 0] - sort_vertexes[0, 0] + 1) //2 < 3 or int(sort_vertexes[0, 0]) > 3:
        if int(sort_vertexes[1, 0]) < 7:
            sort_vertexes[1, 0] = 7
            sort_vertexes[1, 1] = angles[7]

        if int(sort_vertexes[0, 0]) > 3:
            sort_vertexes[0, 0] = 0
            sort_vertexes[0, 1] = angles[0]

    if int(sort_vertexes[3, 0] - sort_vertexes[2, 0] + 1) //2 < 3 or int(sort_vertexes[3, 0]) < 30:
        if int(sort_vertexes[2, 0]) > 28:
            sort_vertexes[2, 0] = 28
            sort_vertexes[2, 1] = angles[28]
        if int(sort_vertexes[3, 0]) < 30:
            sort_vertexes[3, 0] = 33
            sort_vertexes[3, 1] = angles[33]

    # 找出角度差值最大的两条线,即主弯
    new_fix_vertexes = []
    max_angle_vertexes, max_bend_index = max_bend(sort_vertexes)

    if max_bend_index == 1:
        # 主弯是MT, 主弯的下端椎超过腰3就不往下面再找了
        if max_angle_vertexes[1][0] > 27.0:
            # 将主弯MT变成TL/L
            new_fix_vertexes.append(max_angle_vertexes[0, :].reshape((1, 2)))
            tl_angle = round(abs(max_angle_vertexes[0][1] - max_angle_vertexes[1][1]), 2)
            tl_index = int(max_angle_vertexes[1][0])
            for i in range(int(max_angle_vertexes[1, 0]), len(angles)):
                tmp_angle = round(abs(max_angle_vertexes[0][1] - angles[i]), 2)
                if tmp_angle > tl_angle:
                    tl_angle = tmp_angle
                    tl_index = i
            new_fix_vertexes.append(np.array([float(tl_index), angles[tl_index]]).reshape((1, 2)))
            # 从sort_vertexes[0,:] 到 max_angle_vertexes[0, :]搜索PT MT
            down_index = int(max_angle_vertexes[0][0])
            up_kps = kps[0: down_index*2]
            up_angle_index =  calculate_cobb_by_angle(up_kps)
            up_angle_index.sort()
            del up_angle_index[2]
            for i in up_angle_index:
                new_fix_vertexes.append(np.array([float(i), angles[int(i)]]).reshape((1, 2)))

            new_fix_vertexes = np.vstack(new_fix_vertexes)
            new_fix_vertexes = new_fix_vertexes[new_fix_vertexes[:, 0].argsort()]

    # def map_thoracic_index(midle_index):
    #     for k, v in thoracic_index.items():
    #         if v[0] <= midle_index <= v[1]:
    #             return k
    #     return None

    # 生成主弯,代偿弯终板索引的dict
    # compensation_index = 1
    bend_angles = {}

    if len(new_fix_vertexes) > 0:
        cache_vertexes = new_fix_vertexes
    else:
        cache_vertexes = sort_vertexes
    for i in range(len(cache_vertexes)-1):
        # if sort_vertexes[i][0] == extre_peak_valley[0]:
        #     #最大的弯，结构性弯
        #     midle_index = int(np.mean(extre_peak_valley)) // 2 + 1
        #     thoracic_key = map_thoracic_index(midle_index)
        #     bend_angles[thoracic_key] = {'angle': round(abs(sort_vertexes[i][1] - sort_vertexes[i+1][1]), 2),
        #                                 'line_index': [int(sort_vertexes[i][0]), int(sort_vertexes[i+1][0])]}
        # else:
        #     midle_index = int(np.mean(sort_vertexes[:, 0][i: i+2])) // 2 + 1
        #     thoracic_key = map_thoracic_index(midle_index)
        #     bend_angles[thoracic_key] = {'angle': round(abs(sort_vertexes[i][1] - sort_vertexes[i+1][1]), 2),
        #                                 'line_index': [int(sort_vertexes[i][0]), int(sort_vertexes[i+1][0])]}
        #     compensation_index += 1
        thoracic_key = thoracic_keys[i]
        bend_angles[thoracic_key] = {'angle': round((cache_vertexes[i][1] - cache_vertexes[i+1][1]), 2),
                                    'line_index': [int(cache_vertexes[i][0]), int(cache_vertexes[i+1][0])]}

    return angles, peaks, valleys, bend_angles


def test_spine_full_length_photography_noskull_wave():
    import json
    from dataset import sort_points
    from cv_utils import image_show
    root_dir = Path("/home/blake/data/medical/datasets/vertebral/ais_wave")
    phase = 'training'
    save_dir = 'show'

    if not root_dir.joinpath(save_dir).exists():
        root_dir.joinpath(save_dir).mkdir()

    image_items = sorted(root_dir.joinpath('data', phase).glob("*.jpg"))
    for image_file in image_items:
        #wzy.2017.8.28_1.jpg
        #image_file = image_file.with_name("chenxi1_1.jpg") #czh.2018.11.9_1.jpg
        print(image_file.name)
        show_org_image = cv2.imread(str(image_file))
        show_pred_image = np.copy(show_org_image)
        kp_json = image_file.with_suffix(".json")
        angle_json = root_dir.joinpath('data', phase, image_file.stem + '_cobb.json')

        #draw gt kp
        with open(str(kp_json), 'r') as f:
            data_dict: dict = json.load(f)
            landmark = data_dict['shapes']
            landmark, _ = sort_points(landmark, 68)
            # draw kp
            # for i, pt in enumerate(landmark):
            #     show_org_image = cv2.circle(show_org_image, (int(pt[0]), int(pt[1])), 2, (0, 255, 0), 2)
            #     show_pred_image = cv2.circle(show_pred_image, (int(pt[0]), int(pt[1])), 2, (0, 255, 0), 2)

            # for j in range(int(len(landmark) / 4)):
            #     middle_left = ((landmark[j * 4 + 2] + landmark[j * 4]) / 2).astype(int)
            #     middle_right = ((landmark[j * 4 + 3] + landmark[j * 4 + 1]) / 2).astype(int)
            #     cv2.circle(show_org_image, tuple(middle_left), 1, (0, 0, 255), 4)
            #     cv2.circle(show_org_image, tuple(middle_right), 1, (0, 0, 255), 4)
            #     cv2.line(show_org_image, (middle_left[0], middle_left[1]), (middle_right[0], middle_right[1]), (0, 255, 255), 4)
            for j in range(int(len(landmark) / 2)):
                middle_left = (landmark[j *2 ]).astype(int)
                middle_right = (landmark[j * 2 + 1]).astype(int)
                cv2.circle(show_org_image, tuple(middle_left), 1, (0, 0, 255), 4)
                cv2.circle(show_org_image, tuple(middle_right), 1, (0, 0, 255), 4)
                cv2.line(show_org_image, (middle_left[0], middle_left[1]), (middle_right[0], middle_right[1]), (0, 255, 255), 2)

        #draw gt cobb
        # with open(str(angle_json), 'r') as f:
        #     cobb_dict: dict = json.load(f)
        #     angles = [cobb_dict.get("pt", 0), cobb_dict.get("mt", 0), cobb_dict.get("tl", 0)]
        #     lines = cobb_dict.get('lines', [])
        #     # draw selected lines
        #     for line in lines:
        #         cv2.line(show_org_image, (line[0], line[1]), (line[2], line[3]), (255, 0, 0), 6, cv2.LINE_AA)
        #
        #     # draw cobb angles
        #     show_cobb = "GT PT: " + str(format(angles[0], '.2f')) + ", MT: " + str(format(angles[1], '.2f')) + \
        #                 ", TL: " + str(format(angles[2], '.2f'))
        #     show_org_image = cv2.putText(show_org_image, show_cobb, (50, 30), cv2.FONT_HERSHEY_SIMPLEX, 0.5,
        #                                  (0, 255, 0), 1)

        angles, peaks, valleys, bend_angles = calculate_cobb(landmark)


        # draw cobb lines
        for j in bend_angles['main_bend']['vertex_index']:
            middle_left = (landmark[j * 2]).astype(int)
            middle_right = (landmark[j * 2 + 1]).astype(int)
            cv2.circle(show_org_image, tuple(middle_left), 1, (0, 0, 255), 4)
            cv2.circle(show_org_image, tuple(middle_right), 1, (0, 0, 255), 4)
            cv2.line(show_org_image, (middle_left[0], middle_left[1]), (middle_right[0], middle_right[1]),
                     (0, 0, 255), 3)
        for k, v in bend_angles.items():
            if k != "main_bend":
                indexs = v['vertex_index']
                for id in indexs:
                    if id not in bend_angles['main_bend']['vertex_index']:
                        middle_left = (landmark[id * 2]).astype(int)
                        middle_right = (landmark[id * 2 + 1]).astype(int)
                        cv2.circle(show_org_image, tuple(middle_left), 1, (0, 0, 255), 4)
                        cv2.circle(show_org_image, tuple(middle_right), 1, (0, 0, 255), 4)
                        cv2.line(show_org_image, (middle_left[0], middle_left[1]), (middle_right[0], middle_right[1]),
                                 (255, 0, 0), 3)

        show_wave = plot_angles_bar(angles, peaks, valleys)
        show_wave = cv2.cvtColor(show_wave, cv2.COLOR_BGRA2BGR)
        #show_org_image = cv2.resize(show_org_image, (show_wave.shape[1], show_wave.shape[0]))
        show_org_image = mmcv.imrescale(show_org_image, (show_wave.shape[0], show_wave.shape[1]))
        top = (show_wave.shape[0] - show_org_image.shape[0]) //2
        bottom = (show_wave.shape[0] - show_org_image.shape[0]) // 2
        left = (show_wave.shape[1] - show_org_image.shape[1]) //2
        right = (show_wave.shape[1] - show_org_image.shape[1]) // 2
        show_org_image = cv2.copyMakeBorder(show_org_image, top, bottom, left, right, cv2.BORDER_CONSTANT, value=(0, 0, 0))
        show_wave = np.hstack((show_org_image, show_wave))

        put_txt_i = 1
        for k, v in bend_angles.items():
            bend_angles_desc = str(k) + '-->' + 'angle:'+str(v['angle']) + ',' +'index: '+str(v['vertex_index'])
            show_wave = cv2.putText(show_wave, bend_angles_desc, (5, show_wave.shape[0]-put_txt_i*50), cv2.FONT_HERSHEY_SIMPLEX, 1,
                                    (0, 0, 255), 2)
            put_txt_i += 1

        image_show('wave', show_wave)
        cv2.waitKey(0)

        cv2.imwrite(str(root_dir.joinpath(save_dir, image_file.name)), show_wave)


def main():
    print("welcom cobb angles measue!")
    test_spine_full_length_photography_noskull_wave()
    #test_spine_full_length_photography_noskull_pt_mt_tl()


if __name__ == "__main__":
    main()